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Abstract 



The two-basis method to solve the HFB for deformed nuclei in coordinate space 
is examined concerning the precision of the density tail. Small cutoff energies are 
shown to give rise to ripples in the tail, whose wave length corresponds to the cutoff 
■ momentum. More precise solutions require higher cutoff energies, which result in 

HFB matrices of larger dimensions. To circumvent the difficulty of the large dimen- 
sion, we employ another method to solve the HFB — the natural-orbital method 
introduced originally for spherical nuclei — to apply it to deformed nuclei. The 
method can be successfully implemented with a three-dimensional Cartesian mesh 
representation. 



1 Introduction 



In neutron-rich nuclei, the effects of continuum states on the pairing correlation are ex- 
pected to play an important role. As for nuclei near the neutron drip line, it is obvious 
that the continuum states are strongly involved in the pairing correlation because the 
Fermi level A n and the pairing gap A n have similar sizes ~ ±1 MeV. The effect of such 
a continuum-state pairing may be so strong that the neutron drip line can be pushed 
outward by several nuclei. 

For nuclei with A n ~ —5 MeV, too, the continuum states should be taken into account 
explicitly because the pairing-active space for HFB calculations should be larger than one 
major shell, i.e., e cu t > A + hu, in order to take into account a correct size of shell effects. 

The pairing correlation has been treated usually in the BCS approximation, which 
relies on an assumption that the pair-scattering matrix elements ( , 0j'0i|^p|V'iV'j) do not 
depend on the form of the wavefunctions ipi and ipj, e.g. a constant. This assumption 
results in a situation that the nucleus is surrounded by unphysical dilute neutron gas when 
e cu t > 0. Therefore, to include the positive energy states in the pairing correlation, one 
needs to switch from the HF(Hartree-Fock)+BCS to the Hartree-Fock-Bogoliubov (HFB) 
scheme in coordinate space. 

So far, several methods have been presented to solve the HFB equation. For spherical 
nuclei, there are three methods: 

1) Radial differential equations (Dobaczewski et a!., 1984), [0 



2) Finite-element method for finite-range pairing forces (Ring et al., 1997), 

3) Natural-orbital representation (Reinhard et al, 1997). J| 

For deformed nuclei, two methods have been presented. Although spherical cases can 
be solved easily with present computers, deformed cases in coordinate space are still a 
challenge. 

4) Diagonalization in the harmonic-oscillator basis (Gogny et al., 1980), [[| 

5) Two-basis method (Heenen et al., 1994). [|, [7|, || 

In this paper, we describe the formulations in the subsequent two sections and then 
study two subjects: First, in section [|, we examine the two-basis method concerning 
the precision of the low-density tail at large radius as a function of the cutoff energy, 
i.e., the number of HF orbitals included in the diagonalization. We show that more 
precise solutions require higher cutoff energies, which result in HFB matrices of larger 
dimensions. Second, in section |5], we test the natural-orbital method for deformed cases 
in a three-dimensional (3D) Cartesian mesh representation. 



2 HFB in coordinate space 

In this section, we formulate the HF and the HFB in the coordinate-space representation 
in order to elucidate a difficulty of the HFB and to suggest its possible solution in terms of 
the natural-orbital representation. For the sake of simplicity, in this section, we consider 
only one kind of nucleons and designate the number of nucleons by N. The spin of a 
nucleon is represented by s. 

2.1 HF 

In the HF, one should minimize (ty\H\^) for single Slater-determinant states, 

N 

m = IKio>, a) 

S 

by varying {'0i}i=i,-,jV under orthonormality conditions (ipi\ipj) = <%• 



2.2 HFB 

In the HFB, the state takes the following form, 

#basis 

i*> = n w, (3) 

s J 

where "#basis" is the number of basis states of the employed representation: For a 3D- 
mesh representation, |9|] it is the number of the mesh points (times four when spin-orbit 
potentials are included) and typically 10 4 -10 5 . One should vary {(pi, <^i}i=i,... ) #basis under 
appropriate orthonormality conditions. 



The essential difference between the HF and the HFB is that one has to consider 
only iV ~ 10 2 wavefunctions in the former while one has to treat explicitly as many 
single-particle wavefunctions as the number of the basis in the latter. 

2.3 HFB with the two-basis method 

In this method, the HFB equation is solved by diagonalizing the HFB matrix in a single- 
particle basis {ipi}i=i,...,K consisting of the eigenstates of the mean-field hamiltonian h 
(excluding the pairing potential): hipi = e^. The number of the basis, K, is determined 
by a cutoff energy e cut as €\ < ■ ■ ■ < e# < e cu t < £k+i < • • •• We will show in section |] 
that K ^> N to obtain high-precision density tails. 



Figure [I] 



Figure 1: Schematic picture to explain the different number of single-particle states nec- 
essary to express the HFB ground state between the two-basis and the natural-orbital 
methods. 



2.4 HFB in natural-orbital basis 

Owing to the Bloch-Messiah theorem, the state (Q) can be expressed as, 

#basis 

I*) = II (ui + v ia t4) |0>, (5) 
i=i 

a\ = ^2 / dr Tpi(r, s) aj? s : natural orbital (canonical basis). (6) 
s J 

One should vary {ipi,u i }i = i t ...^ asis under constraints for orthonormality, 

{^ j )=6 ij (l<z<j<#basis), (7) 

and for the expectation value of the number of nucleons, 2 v 1 — N. 

Reinhard et al. regard the advantage of the representation (|S|) over (^) to be that one 
has to consider only a single set of wavefunctions {tpi} unlike a double set 
However, we expect more benefit from the natural-orbital representation. Namely, i may 
be truncated as i < K = 0(N) <^ #basis to a very good approximation. It is because ipi 
should be a localized function while the orthogonality does not allow many wavefunctions 
to exist in the vicinity of the nucleus. This situation is illustrated in Fig. |l|. For 3D-mesh 
representations, # basis is proportional to the volume of the cavity (box) while K is 
proportional to the volume of the nucleus. The latter is 10 1 -10 2 times as small as the 
former. 



3 HFB with density-dependent zero-range forces 



3.1 Interaction 



We employ density-dependent zero-range interactions in the rest of this paper for the 
sake of simplicity. There will not be essential differences in the formulation if we use 
the full-form Skyrme force. The force is expressed in the parameterization of the Skyrme 
force as, 

V(n, si; f 2 , Sa ) = (to + \h P S (n - f 2 ) . (8) 

We adopt t = -1099 MeV fm 3 , t 3 = 17624 MeV fm 3+3Q , and a = 0.98 (it is not 1 only 
for the sake of a test of the code) when the force is used to construct the mean-field (HF) 



potential. [ 10| We use a different strength to make the pairing potential. We express the 



force in the parameterization of Ref . Wl 



V p (n, Sl ; f 2 , S2 ) = vj-^- (l - ^p-) 5 (n - f 2 ) . (9) 



We use p c = 0.32 fm (to roughly vanish the volume-changing effect [11 



3.2 Hamiltonian density 

For the sake of simplicity, we treat N=Z nuclei without Coulomb interaction in the rest 
of this paper. Then, the state of the nucleus is expressed as, 



\*) = I[(ui + v i at<4) (ui + Viat<4) |0). (10) 

. , \ / proton \ / neutron 

With the interactions (§|) and @, the total energy for state (|l~0|) is written as, 

E = (m\H\m) = J Ti{r)dr, (11) 

H - ^i^^+H 1 -^ (12) 

where r is the kinetic energy density, 



K 



p(r)=4Y.v*Mr)\ 2 , = 45>^(r)| 2 . (13) 

i=l i=l 

The stationary condition 571 = leads to a mean-filed and a pairing potentials: 

h =f =\v p (l-?-)p. (15) 



P 



c , 



4 Test of the accuracy of the two-basis method for 
HFB 

We have developed from scratch an HFB program for spherical cases based on the two- 
basis method explained in sect. |2.3| . It is used to examine the precision of the low-density 
tail at large radius as a function of the cutoff energy e cut . 



In Fig. ^|we show the density profiles for (Z, N) = (38, 68) calculated with the Braghin- 



Vautherin force ||12|| and v p = —400 MeV. One can see that (i) The density is localized 
with accuracy ~ 10~ 6 fm -3 . (ii) The cutoff makes ripples in the tail. The wave-length of 
the ripples agrees with 27rh(2me cut )~ 1 ^ 2 , i.e., half of the de Broglie wavelength for e cut . 



Figure 2 



Figure 2: Density profiles of HF and two-basis HFB solutions. 

To obtain a more precise density tail, one has to increase e cut , which results in an 
increase of K, the number of single-particle wavefunctions to be considered explicitly. 

1 /2 

K increases only slowly as e C u t for spherical case for each angular momentum, while it 

3 /2 

grows rapidly as e^t for deformed case. The rapid grow of the latter case makes the 
two-basis method practically inapplicable to deformed nuclei because the bottle-necks in 
computation of mean-field methods on 3D meshes are the parts which spend computation 
time proportional to K 2 like the orthogonalization. Therefore, one need a different method 
to solve the HFB for deformed nuclei with enough high cutoff energies. 

5 Implementation of the natural-orbital HFB method 
on 3D mesh 



The natural-orbital method explained in sect. 12^1 was originally introduced for spheri 



cal nuclei. We have implemented the method to treat deformed nuclei in a 3D-mesh 
representation. 

First, let us present a summary of the formulation, which has some differences from 
Ref. H Instead of minimizing E = (^/\H\^/) with given by Eq. (|10|) under constraints 
of Eq. (|7|), one may introduced a Routhian R, 

K K K 

R — E — e Fcrmi • 4 £ v 2 - £ J2 X H - Sij} , Ay = A*,, (16) 

i=l i=l j=l 

and minimize it without constraints. In the above definition, in order to make R real for 
the sake of convenience, K 2 Lagrange multipliers Ay obeying hermiticity are introduced 
instead of \K(K + 1) independent multipliers. Note that 5y is subtracted from (ijj^ijjj) 
in order to treat Ay not as constants like epermi but as functionals of the wavefunctions. 
Stationary conditions of R result in the following equations. 

dR 2 1 1 ha — A ~ 

— = =>- Vi = - - - . assuming ha < 0, (17) 

dVi 2 2 \/(^ - A) 2 + hi 

|£ = o =► Wi-EA w v i -f;f;^{^ i iv*>-^} = o, (is) 

°Vi j=l j=lfc=l 0l ri 

Hi = v 2 h + UiVih. (19) 

For HF, the orthogonalization conditions are easily realized because ipi are eigenstates 
of the same hermite operator h and thus are orthogonal at the solution: The orthogo- 
nalization procedure is needed only because the orthogonality is unstable. On the other 



hand, for the natural-orbital HFB method, the orthogonalization is essential because the 
single-particle hamiltonians TCi differs from state to state. Therefore, the determination 
of the explicit functional form of \j is the most important part of the method. Reinhard 
et al. have proposed 

A y = l^jK^ + n,)^). (20) 



We can justify their choice as follows: From the requirement that Eq. ( |T8| ) must hold at 
the solution (where {ipi\4>j) = <%), one deduces, 

Kj — at the solution. (21) 

Eqs. (^) and fl21~|) are equivalent at the solution because A^- is defined to be hermite. 
Since this hermiticity must hold at any points, one should not adopt Eq. (|2l|) but Eq. (|20[) . 
Because what is needed is Eq. ([21]) and the hermiticity, one may use more complex forms 
like, 

Kj = ~<^-i {Hi + Hj) m {i + h Wi> - (^mA)f) , m 

where f 2 is a parameter to maximize the convergence speed. 

To obtain the HFB solutions in the natural-orbital formalism, one can utilize the 
gradient method, which includes the imaginary-time evolution method: 

5R 

A^A- At— (23) 

We have developed a 3D-mesh natural-orbital HFB program from scratch according 
to the above formulation. First of all, we test the feasibility of the method for 40 Ca. The 
wavefunctions are expressed with 39 x 39 x 39 mesh points with mesh spacing of 0.8 fm. 
Note that the requirement of precision is higher for HFB than for HF because one has 
to treat larger momentum components than the Fermi momentum in HFB. We employed 
the 17-point finite-difference approximation to the Laplacian. The vanishing boundary 
conditions are imposed on the boundary (0th and 40th mesh points) and the wavefunctions 
are anti-symmetrically reflected in the boundary to apply the finite-difference formula. We 
considered K = 20 natural orbitals, which can contain 80 (=2 x A) nucleons. 



Figure 3 



Figure 3: Convergence of HF and HFB in the natural-orbital method. 

We show an example of the convergence history in Fig. |3|. In this calculation, we set 
f 2 = in Eq. fl22|) and At = 10" 24 sec in Eq. ([23]) . We neglect 5\/5ip* in Eq. flT8|). 
Instead, at every 50 steps, {ipi} are Gram-Schmidt orthogonalized in the ascending order 
of ha and then the HFB hamiltonian is diagonalized in the basis to renew {ipi,Vi}. 

In the left-hand portion, the error of Eq. (|18|), i.e., maxj = i ] ... ] ^|Hi'0i — J2j Kjipjl are 
plotted versus the evolution step. The corresponding quantity for HF, maxj =li ... ^/4|^^j — 
{ipi\h\%lji)ipi\ , is also plotted. The figure demonstrates that one can indeed obtain HFB so- 
lutions with the natural-orbital HFB method in the 3D-mesh representation. We obtained 



similar convergence curves for the error of the orthogonality and for the inconsistency be- 
tween the potential and the densities. The right-hand portion shows the error of the total 
energy (estimated as the difference of the total energy from the convergent value), which 
also shows an exponential convergence pattern. 

The speed of the convergence is, however, about ten times as slow as the HF case. We 
are now trying to improve the convergence speed, the stability of the evolution, and the 
robustness of the method. 



6 Influence of pairing correlation on deformation 

As examples of the applications of our natural-orbital HFB program, we show in Fig. [| 
the change of shape due to pairing correlation. The root-mean-square values of x, y, and 
z are plotted as functions of the strength of the pairing interaction v p for 32 S and 60 Zn. 
Both nuclei are triaxial when the pairing correlation is off (pairing gap A = 0). However, 
as soon as the pairing correlation sets in, the axial symmetry is restored. 

The pairing correlation does not always restore symmetric shapes but sometimes break 
a symmetry present without pairing: For 60 Zn, the expectation value of |(1 — t z )Yso is 
zero before the pairing sets in while it is about half the Weiskopf unit when the pairing 
correlation is present. 



Figure 4 



Figure 4: Triaxiality of HFB solutions versus the strength of the pairing interaction. 
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